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Abstract 

We analyze the stability of stationary solutions of a singular Vlasov type hydro- 
dynamic equation (HE) . This equation was derived (under suitable assumptions) as 
the hydrodynamical scaling limit of the Hamiltonian evolution of a system consist- 
ing of a massive piston immersed in an ideal gas of point particles in a box. We find 
explicit criteria for global stability as well as a class of solutions which are linearly 
unstable for a dense set of parameter values. We present evidence (but no proof) 
that when the mechanical system has initial conditions "close" to stationary stable 
solutions of the HE then it stays close to these solutions for a time which is long 
compared to that for which the equations have been derived. On the other hand if 
the initial state of the particle system is close to an unstable stationary solutions of 
the HE the mechanical motion follows for an extended time a perturbed solution of 
that equation: we find such approximate periodic solutions that are linearly stable. 

1 Introduction 

The time evolution of a system consisting of a piston of mass M moving parallel to 
the x-axis in a cube containing non-interacting point particles of unit mass has been 
studied extensively [CLS1, CLS2, CL, G, GF, GP, H, KBM, LPS, Li, PG]. After some 
rescaling of space and time (by the length of the cube) the problem reduces to that of a 
one dimensional system with {Nr) particles in the interval [0, A] (resp., [A, 1]) where 
X(t) is the position of the piston. The left (right) particles move freely between collisions 
with the wall at x = (x = 1) and the piston at x = X(t). At the walls the velocities 
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of the particles get reversed while at x — X(t) the outgoing velocity v' is related to the 
incoming velocity v by the rules of elastic collisions [CLS2, CL, LPS], 



M- 1 



2M 



V 



(1.1) 



v = 



M+l 



v + 



M+l 



where V is the incoming velocity of the piston. It follows from (1.1) that Nl, Nr, as well 
as the total kinetic energy, \ J2f=i vf + \MV 2 (N = N L + N R ) are conserved quantities. 
The dynamics of the system can be reduced to a billiard in a (2N + l)-dimensional 
domain (polyhedron), cf. [CL]. It was shown in [LPS, CLS1, CLS2], under certain quite 
restrictive conditions on the initial distribution of gas particles, that, in the limit N — > oo, 
M ~ iV 2 / 3 , the dynamics of the piston and the gas satisfy a closed system of Euler-type 
hydrodynamical equations (HE) for a time interval (0, r) in which any particle had at 
most two collisions with the piston. 

The origin of the scaling M ~ iV 2 / 3 is as follows. For N particles with velocities of 
0(1) distributed with density of 0(1) in a parallelopiped of length L and crossectional 
area A the number of particles colliding with the piston per unit (unrescaled) time and 
hence the pressure (from each side) is proportional to A. To ensure that, on this time 
scale, the acceleration of the piston stays of 0(1) as L, A and iV ~ O(AL) grows to 
infinity it is necessary to make the mass of the piston grow as A. For a cube this 
corresponds to M ~ iV 2 / 3 . In the rescaled units the number of collisions experienced by 
the piston per unit time is of 0(N) independent of A and the HE hold for general M as 
long as M ~ N a , a G (0, 1), i.e. when the kinetic energy of the piston becomes negligible 
compared to that of the gas. The time interval (in the scaled units) during which the 
derivation of the HE remains valid depends on a (getting larger as a — > 1), see Remarks 
3 and 4 in Sect. 4 of [CLS2]. For a = 2/3 this time interval is such that the piston 
suffers no more than two collisions with any gas particle. It is however not clear from the 
derivation to what extent those equations may actually approximate the real evolution 
of the system with large N, for longer times. This led us to carry out extensive computer 
simulations of particle systems, with M ~ iV 2//3 , iV as large as 27 x 10 6 [CL], and initial 
conditions for which the hydrodynamic equations have trivial solutions X(t) = 0.5 and 
V(t) = for all t > 0. We found nevertheless that for certain initial velocity distributions 
the trajectory of the piston diverged from these values. In particular it was observed in 
these simulations that the particle plus piston system, after experiencing large random 
fluctuations, quickly converges to a more stable regime, in which the piston and the gas 
undergo regular slowly damped oscillations lasting a long time. The parameters of those 
oscillations (the period, the amplitude, and the rate of damping) seem to depend little 
on many details of the initial distributions. 

A possible explanation for this behavior is that the mechanical trajectory of the system 
with iV large but finite, being subjected to the intrinsic random fluctuations associated 
with the discrete nature of the gas particles, behaves as a perturbed solution of the HE, 
i.e. as if the initial conditions were subjected to a small random perturbation. There 
large fluctuations in the particle dynamics observed experimentally in [CL] may thus be 
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due to the instability of the HE for the initial conditions considered. This would imply 
that, on the contrary, when the solution of the hydrodynamical equations is stable, then 
the mechanical evolution of the system should remain close to that solution for a very 
long time. 

Here we present evidence, but no proof, for this conjecture. We first find a class of 
initial conditions for which the solution of the HE are stable and note that simulations 
with such distributions in [CL] indeed yielded a piston trajectory, which remained close 
to the solution of the HE with V(t) ~ . Furthermore, we observe that the mechanical 
piston trajectories, which do not stay close to the solution of the HE with the prescribed 
initial conditions, appear to follow some perturbed oscillating solutions of the HE equa- 
tions. These solutions appear to act as an "attractor" for the HE 1 . We cannot rigorously 
prove the existence of attracting periodic solutions for the HE yet, but we present some 
approximate constructions of such solutions. 

The paper is organized as follows. In Section 2 we state the hydrodynamical equations. 
In Section 3 we prove global stability for a class of stationary solution. In Section 4, we 
use a perturbative analysis to find sufficient conditions for linear instability. In Section 
5 we investigate a particular family of stationary solutions which include those used in 
the [CL] simulations and show that it contains both linearly stable and unstable ones, 
alternating in a very intricate manner. In Section 6, an approximate construction of 
periodic solutions of the HE is presented. 

2 Hydrodynamical equations 

Let X(t) G (0, 1) be the position of the piston at time t and V(t) its velocity. We denote 
the continuum density of the gas in [0, 1] x R by a function p(x, v, t). The HE describing 
the time evolution of this, continuum fluid plus piston system, are as follows. 

(HI) Free motion. Inside the container the density satisfies the standard continuity 
equation for a noninteracting particle system without external forces: 

for all x G (0, 1) except x — 0, x — 1 and x = X(t). 
Equation (2.1) has a simple solution 

p(x, v, t) — p(x — vs, v, t — s) (2.2) 

for < s < t such that x — vr^ {0, X(t — r), 1} for all r G (0, s). Equation (2.2) has one 
advantage over (2.1): it applies to all points (x,v), including those where the function p 
is not differentiable, or even continuous. 

1 Duc to an obvious time-reversibility of the hydrodynamical equations, there can be no attractors, 
strictly speaking, but saddle-type periodic solutions can very well exist. 
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(H2) Collisions with the walls. At the walls x = and x = 1 we have 



p(0,v,t)=p(0,-v,t) (2.3) 

p(l,v,t)=p(l,-v,t) (2.4) 

(H3) Collisions with the piston. At the piston x — X(t) we have (this is obtained from 
(1.1) when M -> oo) 

p(X(i)-0,u,t) = p(X(*) -0,2V(t) for u<V(t) 

p(X(*) + 0,u,f) = p(X(t) + 0,2V(t) -v,t) for v>V(t) (2.5) 

where v represents the velocity after the collision and 2V(t) — v that before the 
collision; and 

X(t) = X(0) + / V(s)ds (2.6) 
Jo 

is the (deterministic) position of the piston. 

It remains to describe the evolution of V(t) (which we take to be left continuous). 
Suppose the piston's position at time t is X and its velocity V. The piston is affected 
by the fluid at (x, v) exerting pressure on it from the right (x = X + and v < V) and 
from the left (x = X — and v > V). Accordingly, we define the density of the fluid in 
contact with the piston ( "density on the piston" ) by 

q (v,t;X,V) = \ P % + °> V >% ' l ' V< l (2.7) 
^ v ' ' ' ' \ p(X - 0, v, t) if v > V v ' 

(H4) Piston's velocity. The velocity V = V(t) of the piston satisfies the equation 

{v - Vf sgn(w - V) q(v, t; X, V) dv = (2.8) 

— oo 



/oo 
-oo 



The origin of eqs. (H1)-(H3) in the particle system is clear. Eq. (H4) is essentially a force 
balance equation — since the rate of collision of the piston with particles on either side and 
consequent force on piston is much larger than mass of the piston when N/M — > oo, V 
adjusts instantaneously to make the forces from the two sides balance exactly. The system 
of (hydrodynamical) equations (H1)-(H4) is now closed and, given initial conditions, 
satisfying (2.3)-(2.8) at t — 0, completely determine the functions X(t), V(t) &ndp(x, v, t) 
for t > 0. When the initial conditions do not satisfy these equations one has to imagine 
that they become satisfied instantaneously for t — 0+. The existence and uniqueness of 
solutions of (H1)-(H4) were proven, under general conditions, in [CLS1, CLS2]. We need 
only to assume that the p(x, v, 0) is bounded, piecewise differentiable, and either has a 
compact support in the x, v plane or decays fast enough as \v\ — > oo. We also require 
that / p(x, v, 0)dv > for all x. 
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The HE, like the Vlasov equations for plasmas, are time-reversible, see [P] and [MP] . 
They preserve the classical integrals of motion. The mass of the fluid to the left and to 
the right of the piston as well as the total kinetic energy of the fluid remain constant 
along any solution. 

(Dl) Mass conservation 

Ml — [ [ p(x,v,t) dv dx, Mr= [ [ p(x,v,t) dvdx 
Jo J Jx(t)J 

(D2) Energy conservation 

2E= f [ v 2 p(x,v,t)dvdx (2.9) 



Equation (2.8) also preserves the total momentum // v p(x,v,t) dv dx, but the latter 
changes due to reflections at the walls. We note that the piston itself does not contribute 
to the total momentum and energy of the system in this model, because its mass and 
energy vanish, when divided by N, in the limit N — > oo. (The mass, energy and momen- 
tum of the fluid all correspond to the original quantities in the particle system divided 
by TV). 

The HE define a dynamics on the domain G := {(x,v) : < x < 1} in which every 
point (x, v) G G moves freely with constant velocity and collides elastically with the walls 
and the piston. Denote by (x t ,v t ) the position and velocity of an arbitrary point at time 
t > 0. Then (HI) translates into x t = v t and v t — whenever x t ^ {0,1, X(t)}, (H2) 
becomes (x t +o,v t +o) = (xt-o, —Vt-o) whenever x t -o G {0, 1}, and (H3) gives 

(x t+0 ,v t+ o) = (x t - ,2V(t) - vt-o) 

whenever x t -o = X(t). Note that the point (x t ,v t ) moves in G and reflects at the walls 
and the piston as if those had infinite mass. 

The motion of points in G is described by a one-parameter family of transformations 
F* : G — > G defined by F t (xo,Vo) = (x t ,v t ) for t > 0. We will also write F~ t (x t ,v t ) = 
(x ,v ). According to (H1)-(H3), the density p(x,v,t) satisfies a simple equation 

p(x t , v t , t) = p{F~\x t , v t ), 0) = p(x , v , 0) 

for all t > 0. It is easy to see that for each t > the map F t is one-to-one and preserves 
area, i.e. det \DF t (x,v)\ = 1. Hence, the family F l describes an incompressible flow on 
G and consequently: 

(D3) Incompressibility. For any a < b the Lebesgue measures (areas) of the sets 

{(x, v) : a < p(x, v,t) < b, < x < X(t)} 

and 

{(x, v) : a < p(x, v, t) < b, X(t) < x < 1} 
remain constant in time 
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A particular case in which it is possible to solve equations (H1)-(H4) analytically, 
is when the initial distribution is stationary. This happens when p(x, v,0) satisfies two 
conditions: 



(SI) Uniformity and symmetry. The initial density p(x, v, 0) = p(x,v) is of the form 

p(x,v) = 



Pl(\v\) for x < X 
Pr(\v\) for x > X 



for all v and X(0) = X . 
(S2) Pressure balance. The pressure on the piston from both sides is equal: 

P L := 2 / v 2 p L (v) dv = P R :=2 v 2 p R {v) dv (2.10) 

JO JO 

Under conditions (S1)-(S2) the equations (H1)-(H4) have a simple solution: the 
system remains frozen in its initial state: 

X(t)=X , V(t) = 0, p(x,v,t)=p(x,v,0) (2.11) 

for all t > 0. 

We will analyze in the next three sections the stability of this stationary solution. 
Note that there is no requirement on the form of Pz,(H) or Pr(\v\); all that is required is 
a balance of forces (2.10). 



3 Globally stable solutions 

Here we consider stationary solutions p(x, v) satisfying (S1)-(S2) and an additional mono- 
tonicity requirement: 

Pl(N) > Pl(\v 2 \) and p R (\ Vl \) > p R (\v 2 \) (3.1) 

for all \vi\ < \v 2 \. We claim that such solutions are globally stable. This criteria is very 
similar to the stability criteria for the Vlasov equation described by Penrose [P] and by 
Marchioro and Pulvirenti [MP]. 

Before we state our result, we introduce some notation. Denote by || ■ || the following 
special norm on the space of functions on G 

\\f(x,v)-g(x,v)\\ =JI \f(x,v)-g(x,v)\(l + v 2 )dvdx (3.2) 

Theorem 3.1 Let p(x,v) satisfy (SI), (S2) and (3.1). Then for any e > there exists 
5 > such that if the initial density p(x,v,0) satisfies \\p(x,v,0) — p(x, v)\\ < 5, and 
X(0) = X then 
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(i) \\p(x,v,t)-p(x,v)\\ < e; 



(ii) \X(t)-X(0)\ < e 
for all t > 0. 

Proof. The proof of claim (i) is identical to the proof of the stability theorem of Marchioro 
and Pulvirenti [MP]. We note that their theorem is stated in the L 1 norm, but, in fact 
it was proven in the (3.2) norm. (There is a minor error in [MP], the conditions that 
assure that the Li norm and the (3.2) norm are equivalent are not explicitely assumed.) 

It is clear from (2.9) that, given the position of the piston X and values of the areas 
of the level sets, defined in (D3), the minimal possible value of the total energy for any 
phase-space density ir(x,v), is attained when ir(x,v) is uniform in x and monotonically 
decreasing in \v\ in each compartment. 

Consider first the case where tt(x, v) has the same area of the level sets as some p(x, v) 
satisfying (SI) and (S2). Then the minimum of the energy when the piston position is 
X is attained when 

tt(x,v) =p L (vX/X Q ), 0<x<X 

and 

n(x,v)=p R (v(l-X)/(l-X )), X<x<l 
The minimal total energy is then 

roo rX roo rl 

2E min (t) = v 2 ir(x,v) dxdv + / / v 2 ir(x, v ) dxdv 

Jo Jo Jo J x 

(we used a change of variable u = vX/X in the first integral and u = v(l — X)/(l — X ) 
in the second one). Using the pressure balance (2.10) and denoting P = P L = P R gives 

E (t)- P ( X ° I t 1 -*^ 

min[ ) 2 \xi + (i-xy ) 

Consider now the above expression as a function of X. Its minimum is attained at the 
point where dE m { n /dX = 0, i.e. 

X 3 (1-X) 3 

which is only possible if X = X . Therefore, the state X = X provides a unique 
minimum of the total energy function under the incompressibility constraint (D3). This 
proves the theorem when the p(x,v, 0) has exactly the same area of the level sets as p(x, v). 
For perturbed initial densities p(x,v,0) the above estimates only hold approximately, 
hence small changes in the system are possible, but large changes are not. This proves 
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the theorem. Slightly extending the above theorem, let us assume that the initial density 

p(x, v,0) satisfies (SI), (3.1), but not (S2). If the pressures Pl and Pr in (2.10) differ 
by a small amount A = \Pl — Pr\, then one can estimate how far the piston can swing. 
The piston can move as long as 



By simple calculations one obtains, to the leading order of A, the following bound on 
the piston displacements: 



4 Perturbative analysis 

Here we consider initial densities p(x, v,0) satisfying (S2) and the following stricter ver- 
sion of (SI): 

(S3) Full uniformity and symmetry. The initial density p(x, v,0) is uniform in x across 
the entire cylinder, i.e. p(x,v, 0) = £>o(M) for all v and < x < 1. 

We also assume that the piston is initially at the midpoint X(0) = 0.5. Of course, 
under the conditions (S2)-(S3), the hydrodynamical equations (H1)-(H4) have a simple 
stationary solution (2.11). On the other hand, we no longer assume monotonicity (3.1). 
We use perturbative analysis to investigate the stability of the stationary solution (2.11). 

From now on we denote by po(v) = Po(\v\) an initial density satisfying (S2)-(S3) and 
by p(x, v, 0) a perturbed initial density, which we write down as 

p(x, v, 0) =po(\v\) + epi(x, v, 0) 

where epi(x, v, 0) is a small perturbation. We will work to first order in e, i.e. ignore 
terms of order o(e). For t > 0, we decompose the density p(x,v,t) as 

p(x, v, t) = po(\v\) + £Pi(x, v, t) 

We also set p±(x, v, t) = Pl(x, v, t) for x < X(t) and pi(x, v, t) = Pr(x, v, t) for x > X(t). 
According to (2.8), the velocity V(t) of the piston is given by 




v 2 p(X(t) + 0,v,t)dv < 2E(0) 



where X = X(0), as before, and E(0) is the initial total energy: 
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where X = X(t) is the position of the piston. Expanding in e and ignoring terms of 
order o(e) gives 

/ °° v 2 p L (X, v, t) dv - J ^ v 2 p R (X, v, t) dv 
U £ 4J °°vp (v)dv 

Note that by integrating by parts we obtain (for piecewise smooth p ) 

2 / vpo(v)dv = — / v 2 p'Jv) dv 
Jo Jo 

We define 

h(v) = —p' (v) for v > 
and for the sake of completeness, set h(—v) = h(v). Then 

T/m Jo°° v 2 p L (X, y, t) dv - ^(X, ^ t) dv 
V{t) = 6 2f °°v 2 h(v)dv 

The function po(v) does not have to be differentiate, we interpret —h(v) here as the 
generalized derivative of po{v). Denote by (•)+ the integration J °° -dv, and by (•)_ the 
integration J ^ ■ dv. Then 

v( ]= (v 2 p L (X,v,t)) + - (v 2 p R (X,v,t))^ 
[) £ 2{v 2 h{v)) + 

The density of the gas after interaction with the piston is given by the formulas (2.5), 
which imply: 

p(X-0,-v,t) = p(X - 0, v + 2V, t) 

= p (v + 2V)+sp L (X,v,t) 

= p (v) + 2Vp' (v)+ep L (X,v,t) 

= p (v)-2Vh(v)+ep L {X,v,t) 

Here we assume v > and ignore terms of order o(e). Hence, the "reflection rule" can 
be written as 

, Y *\ (V + \ U \ (v 2 PL(X, V,t)} + - (v 2 p R (X,V,t)}_ 

p L {X, -v, t) = p L {X : v, t) - h(v) 



(v 2 h(v)) + 

This expression suggests the introduction of new functions: 

q RtL (x,v,t) = — 

h(v) 

and 

v 2 h{y) 

p(v) = 



yh(v)). 
9 



The above expression for Pl(X, —v,t) can now be written as 



q L {X, -v, t) = q L (X, v, t) - (q L (X, v, t) p{v)) + + (q R (X, v, t) p(v))- (4.2) 

Similarly, on the other side of the piston, 

q R (X, v, t) = q R (X, -v, t) + (q L (X, v, t) p(v)) + - (q R (X, v, t) p(v))- (4.3) 

One can interpret these "reflection rules" as follows: the functions qi and q R "exchange" 
their average values with respect to the "density" p{v). 

Note that p{v ) is normalized, so that (p(v)) + = 1, but it is not necessarily positive (or 
even nonnegative). On the other hand when p{v) > 0, i.e. the unperturbed density po( \v\) 
is nonincreasing, thus satisfying (3.1). In this case the stationary solution (2.11) is stable, 
as we already know by Theorem 3.1. Here we recover this result by our perturbative 
analysis: 

Theorem 4.1 The quantity 

Q = jj q2{x ^ p{v) d X dv 

is constant in time, i.e. dQ/dt = 0. Here q = qi for x < X and q = q R for x > X . 

Proof. Clearly, Q cannot change just due to the free motion of the gas or due to collisions 
with the walls, so we only need to worry about collisions with the piston. The gas particles 
colliding with the piston during an infinitesimal interval (t, t + dt) lie in the two triangles 
on the xv plane: X — v dt < x < X for v > and X < x < X — v dt for v < 0. The 
outgoing particles lie in similar symmetric triangles. Hence, during the interval (t,t + dt), 
the quantity Q decreases by (up to the factor of dt) 

r°° ql(X,v,t)p(v) ^ | ro q R (X,v,t)p(v) ^ 

JO \v\ J-oo \v\ 

= (q 2 L (X, t) p(v))+ + (q 2 R (X, v, t) p{v)). 

and it increases by 

gL{X,v,t)p{v) ^ | r°Q <&(X,v,t)p{v) dv 

J-oo \v\ Jo \v\ 

= (q 2 L (X, v, t) p{v))_ + (q 2 R (X, v, t) p(v))+ 

After substituting (4.2) and (4.3) into the above expressions for and q R and some 
manipulations, all changes in Q cancel out and so it stays constant. □ 

When Po(M) is strictly decreasing, hence p(v) > 0, then Q is a norm in the space of 
functions. Thus, the above theorem implies linear stability. 



10 



When po(M) is decreasing, but not strictly, then p(v) > 0, but there may be regions 
where p(v ) = 0. They correspond to the intervals where p' Q = 0, i.e. where po is constant. 
On such intervals, the reflection rules (4.2)-(4.3) for the perturbed density Pl,r are trivial: 

p L (X, -v, t) = p L (X, v, t) and Pr(X, v, t) = p R (X, -v, t) 

In this case pi and p R cannot grow either. Therefore, we obtain linear stability for all 
nonincreasing po ( | v | ) . 

Next we turn to unstable solutions. The stationary solution for an initial density 
Po(M) satisfying (S2)-(S3) is linearly unstable if some small perturbations grow expo- 
nentially in time, i.e. v, t)\\ ~ A* for some p±(x, v,0) and A > 1. This is equivalent 
to having a positive Lyapunov exponent in the subspace spanned by the function p\ 
and its images. To investigate the existence of such perturbations we first simplify the 
collision rules (4.2) and (4.3). Consider the following "symmetric" and "antisymmetric" 
linear combinations of q^ and q R . 

q+(x, v, t) = [q L (x, v, t) + q R (l - x, -v, f)]/2 

and 

q_(x,v,t) = [q L (x,v,t) - q R (l - x,-v,t)]/2 
They are defined for x < 1/2. The collision rules (4.2)-(4.3) now take form 

q + {X,-v,t)=q + (X,v,t) (4.4) 

and 

q.(X, -v, t) = g_(X, v, t) - 2(g_(X, v, t) p{v)) + (4.5) 

Hence q + is simply a periodic function in t, so it cannot grow to infinity or decrease to 
zero. In other words, it cannot affect the stability or instability of the hydrodynamical 
equations. The latter is determined by g_ alone. So we will only consider q_ and omit 
"— " for brevity. Our collision rule then reduces to a single equation 

q(X, -v, t) = q(X, v, t) - 2(q(X, v, t) p(v)) + (4.6) 

Next we demonstrate, by example, that densities Po(M) f° r which the stationary so- 
lution (2.11) is unstable do exist. 



Example. Let p be a rectangular function defined by 

, v f 1 if 0.5 < \v\ < 1 
P ^ V) = { otherwise (47) 

This po(v) satisfies (S2)-(S3) but not (3.1). We will show that the corresponding sta- 
tionary solution is linearly unstable. 
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First, the function h = 



—p' is the sum of two delta functions: 



h(v) = -<5 . 5 + 5i 
(and symmetrically for v < 0). It is easy to compute p directly 

1 r 4 r 

P = "J ^0.5 + g Ol 

Now the reflection rule (4.6) implies 

?(-l) = -| + | 9(0.5) 
g(-0.5) = -|g(l) + |g(0.5) 

Note that only the values ±0.5, t) and p(x, ±1,£) will evolve in a nontrivial way, as 
specified above, since h(v) = for all v ^ {1, 0.5, —0.5, —1}. 

We now construct a linear subspace of functions q = g_ that stays invariant under 
the above transformations and in which functions grow exponentially in time (since the 
q + component of the perturbation is irrelevant, we set it to zero). We can simplify the 
construction further by assuming that at time t — 

q(x, ±1, 0) = u±, q(x, 0.5, 0) = u 2 , q(x, — 0.5, 0) = u 3 

with some constants 1*1,1*2, 1*3 (the choice of indices 1,2,3 is rather arbitrary). We note 
that the functions Pl,r(x, v, 0) are now piecewise constant and are completely described 
by the values 1*1,1*2,1*3- The space of such perturbations is three-dimensional. 

It is easy to see that at time t — 1 the functions j>l,r will again be constant on the 
same intervals, hence they will be described by some other constants 1*^,1*3,1*3. Our 
collision rule (4.6) implies that the vectors u' = (u[, u' 2 , i* 3 ) T and u = (1*1, 1*2, U3) T are 
related by a linear transformation 

u' = Au 

where A is a 3 x 3 matrix: 

A=i 

3 

After that, the evolution will proceed periodically - the vector u will be multiplied by 
the matrix A at times t — 1, 2, 3, . . .. The matrix A has three real eigenvalues: 

-4±y/7 

Ai, 2 = and A 3 = 1 

The largest eigenvalue A = —(4 + v / 7)/3 ~ —2.215 has the following (unit) eigenvector: 

u = (0.4472, -0.3680, 0.8152) 



/ -5 


2 


\ 








3 
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This eigenvector spans a one-dimensional subspace in the space of perturbation densities, 
which is invariant over the time period t — 1 and in which the corresponding perturba- 
tions are expanded by a factor |A| ~ 2.215. Roughly, the perturbations double over one 
period. 

To explore the above periodic growth of perturbations, we note that the piston velocity 
is given by 

V(t) = £ -((q L (X,v,t)p(v)) + -(q R (X,v,t)p(v))_) 
= e(q(X,v,t)) + 

Hence in our example, during the time interval < t < 1 

V = - (4«i -u 3 ) = 0.9736 e 
During the next time interval 1 < t < 2 we have 

V = E - (Au[-u' 3 ) = -2.156 £ 

o 

and so on. Hence, over a unit period of time, the piston velocity grows by a factor of ]A| = 
2.215 and changes sign - the piston starts its movements back and forth (oscillations) 
that increase exponentially in time. 

We note that the same density (4.7) was studied in [CL] where the trajectory of 
the piston was computed after an initial configuration of gas molecules was selected 
randomly from the distribution po(v) given in (4.7). It was found [CL] that the piston 
indeed made oscillations which increased exponentially in time. The piston's velocity 
grew as const-/?* with R m 1.6. This experimental estimate is in a reasonable agreement 
with our calculation of the largest eigenvalue 2.215. 

Next, we modify the unstable perturbations q found above and make them smooth 
(rather than piecewise constant) functions of v. 

We will be looking for the function q of the form 

q(x,v,t) = C(v)e z{t ' x/v) 

where z is a complex constant. Note that due to (2.2) the function q (with v fixed) 
can only depend on t — x/v. We chose the exponential form in order to investigate the 
existence of solutions of the linear equation which grow exponentially with time. Also, 
for convenience, we introduce the new space coordinate y in the following way: for all 
v > and x < 0.5 we set y = x + 0.5, for v < and x < 0.5 we set y = 0.5 — x, for 
v > and x > 0.5 we set y = x — 0.5, and for v < and x > 0.5 we set y — 1.5 — x. The 
so designed coordinate y assumes value zero when a point (x, v) e G moving under F* 
reflects off the piston, then grows to 0.5 when the point travels to the wall, and grows 
further from 0.5 to 1 when the point travels from the wall back to the piston. 
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In the new coordinate y, we will be looking for the function q of the form 

q(y,v,t) = C{\v\)e z ^- y ^ 

More precisely, let 

q(y,±l,t) = C(l)e?M 

q{y, ±0.5, t) = C{0.5)e z(t ' 2y) 

Recall that po(\v |) is the characteristic function of the interval [0.5, 1]. 
Now, the reflection rule (4.6) implies 

C(l) = -^(l)e- 2 + ^(0.5)e^ 

C(0.5) = -^C(l)e- 2 + ^C(0.5)e-^ 

We need to find z for which the above system of equations has a nontrivial solution. Put 
e z = A and introduce an auxiliary variable -D(0.5) = C(0.5) e~ z . Now the above system 
can be rewritten as 



AC(l) : 

AD(0.5) = C(0.5) 



-|C7(1) + |I>(0.5) 



-_C(l) + -£>(0.5) 



AC(0.5) : 

Hence, A is an eigenvalue of the matrix of coefficients 



-5 2 \ 
3 
-8 5 / 



which is the same matrix A that we encountered before. We take its leading eigenvalue 
|A| > 1 and set 

z = In |A| + in 

The function q now takes form 

q(y, v, t) = ±C(\v\) |A|* _2//H cos n(t - y/\v\) 

where C(0.5) and C(l) are the coordinates of the leading eigenvector, and we only take 
the real part, for obvious reasons. Since |A| > 1, we have an exponential growth of 
perturbations and thus linear instability. This gives us smooth unstable perturbations. 

We now generalize the above construction to arbitrary nonmonotonic initial densities 
p . Let po(v) satisfy (S2)-(S3) but not (3.1). We will be looking for perturbations of the 
form 

q(y,v,t) = C(\v\)e z{t - yM) (4.8) 
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with the same convention on y as before. The reflection rule (4.6) leads to (cancelling 

e zt ) 

Civ) = Civ) e- z/v - 2 f° Civ) e- z/v p(v) dv 

Jo 

for all v > 0. Denoting 

D = -2 J C(v) e- z/v p{v) dv 

gives immediately 

0(.)= D 



1 - e~ z / v 

Thus, we not only eliminated t but determined the function C(v) up to a constant factor. 
The above solution exists if 

D = -2 

Jo 1 - e- z /'° 

or, cancelling D, 

To 1 - e 2 A- 2 v ; 

If this equation has a solution z with Re(z) > 0, we immediately obtain an unstable 
perturbation (4.8). Otherwise our construction of unstable perturbations does not work. 

Unfortunately, it does not seem to be easy to solve equation (4.9) for particular 
functions p(v) or even to determine if it has solutions with a positive real part, as we will 
demonstrate in the next section. 

Next we mention an important property of (4.9). Let us denote 

F(z) := I" - P{V \ dv- 1 - (4.10) 
v ; Jo l-e z l v 2 v ; 

Lemma 4.2 F(z) + F(-z) = for all z G C. 
Proof. 

F(z) +F (-z) = r^dv- r e ^^dv-i 

K J K J Jo l-e z ' v Jo l-e z / v 
r oo (! _ e z / v )p(v) 



r°° (l - 

Jo r 



dv-1 



I) □ 



As a result, the existence of a solution of (4.9) with Re(z) > is equivalent to that of 
a solution with Re(z) < 0. The alternative is when all the solutions lie on the imaginary 
axis Re(z) = 0. 
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5 A special family of densities 

Here we investigate a family of rectangular densities 



Po{v) 



1 if r < \v\ < 1 
otherwise 



(5.1) 



where < r < 1 is the parameter of our family. Note that our example (4.7) is a 
particular case of (5.1) with r = 1/2. It is easy to compute 



and 



p(v) 



h(v) = -p' (v) = 6 r (v) - 6x(v) 
v 2 h(v) 



1 



/ °° v 2 h(v) dv l-i 

Since h(v) = for all v £ {1, r, — 1, — r}, we only consider perturbations q(x, v, t) defined 
for v = l,r, —1, — r. The reflection rule (4.6) now gives 



q{X,-l,t) 
q(X,-r,t) 



-aq(X,l,t) +(3q(X, r,t) 
-jq(X,l,t) + aq(X,r,t) 



where 



1 +r 2 



a 







2r u 



7 



(5.2) 



While we cannot construct unstable perturbations for arbitrary irrational values of r, it 
is relatively easy to investigate the case of rational r. But even in this case, we can only 
provide partial answers leaving some questions open. 



P m-1 



P 3 



1/n 



P m+1 



P 2 1/m Pj 



Figure 1: The construction of points Pi. Here m = 7 and n = 12. 
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Let r = m/n be a rational number with 1 < m < n. Now (5.2) takes form 

2m 2 2n 2 



n 2 + m 2 



a = 



rr 



7 



(5.3) 



To investigate the evolution of perturbations q(x,v,t) as t grows, we consider n + m 
points Pi G G, 1 < i < n + m, shown on Fig. 1. The points Pi are defined as follows: 

(-1,0.5 - (i - l)/m) forl<i<m/2 + l 
(1, (i — l)/m — 0.5) for m/2 + 1 < % < m 

(— r, 0.5 — (i — m — l)/n) for m < i < m + n/2 + 1 



r, u 



m — l)/n — 0.5) for m + n/2 + 1 < « < m + n 



It is crucial to observe that the points Pi move under the dynamics in a periodic fashion. 
In a time period At = l/m, the point Pi is mapped to P i+i for all 1 < i < m and all 
m + l<i<m + n. Also, P m moves to the piston, gets reflected off it and lands on P ± . 
Likewise, P m +n moves to the piston, gets reflected off it and lands on P m+1 . Therefore, 
the time shift At permutes the points Pi, 1 < % < m + n, in two independent cycles. The 
reason why we combine the two cycles together is that they are linked by the reflection 
rule, as we will see shortly. 

For each i, denote %(t) = q(Pi,t). Then we have 

qi(t + At) = q i - 1 (t) 

for all 2 < i < m and m + 2<i<m + n. The reflection rule now implies 

gi(t + At) = -aq m (t) + f3 q m+n {t) 
q m+1 (t + At) = -7 q m (t ) + a q m+n (t ) 

Thus, the vector q(t) = (?i(t), . . . , q n +m(t)) is updated at time t + At by the rule 

q(t + At) = B q(t) 

where B is an (n + m) x {n + m) matrix, 

/0 ••• 



B 



1 



—a 



: : 

1 
-7 

1 



a 



1 



17 



We conclude that the existence of unstable perturbations q(t) is equivalent to the exis- 
tence of an eigenvalue A of B such that |A| > 1. The characteristic polynomial of the 
matrix B is 

P(A) = \ m+n + a\ n -a\ m -l (5.4) 
where a = (n 2 + m 2 )/(n 2 — m 2 ) is defined in (5.3). 

Remark. Interestingly, the equation (4.9) can be reduced to -P(A) = as well. Indeed, it 
is easy to see that 



f 

Jo 



»<"> 1 



^2 



z/r 



Now the substitution A = e z l m and some algebraic manipulations show that Eq. (4.9) is 
equivalent to P(\) = 0. 

It is easy to see that if A is a root of -P(A), then so is 1/A (this reciprocability also 
follows from Lemma 4.2). Thus, the existence of unstable perturbations is equivalent to 
the existence of eigenvalues of B that do not lie on the unit circle |A| = 1. 

If an eigenvalue |A| > 1 of B exists, then the perturbations in the corresponding 
eigenspace grow by the factor of |A| over the time period At = 1/m. Hence, the expansion 
factor per unit time would be A = |A| m . 

Theorem 5.1 Let r = m/n be a rational number with an even denominator n (hence, 
m is odd). Then there is a unique eigenvalue of~B such that A < — 1. This eigenvalue has 
multiplicity one. The expansion factor per unit time A r = \X\ m depends on r continuously, 
and we have, asymptotically, 

A r = 1 + const • r 3/2 + C(r 2 ) as r^O (5.5) 

and 

const 1 . „. 

A r ~ as r — > 1 (5.o) 

1 — r 



Proof. One can easily check that, under the conditions of the theorem, P(— 1) > and 
P{— oo) < 0, hence a root A < — 1 exists. Next, 

P'(A) = [(n + m)\ n + an\ n - m - am] A" 1 " 1 

and so P'(-l) < and P'(-oo) > 0. Now let Q(A) = (n + m)\ n + an\ n ~ m - am, then 

Q>(\) = [( n + m)\ m + a(n - m)] n\ n - m ^ 

and so clearly Q'(\) < for all A < 1. Putting these facts together proves the uniqueness 
and the simplicity of the root A < 0. 



18 



The equation P(X) = can be rewritten in terms of A r = |A| m as follows: 

-K +1/r + \ ± ^ 2 K /r + \ ±J ^ 2 K -1=0 (5.7) 

Now the continuity of A r , as a function of r, is obvious. Note that our argument is only 
valid when r = m/n with an even n and an odd m, because this parity condition dictates 
the signs in (5.7). 

To prove (5.5), one can substitute A r = 1 + e into (5.7) and expand all the terms in 
Taylor series, the calculation is then straightforward and we omit it. The proof of (5.6) 
is similar. □ 

J In A 




0.2 0.4 0.6 0.8 1 



Figure 2: logA r as a function of r. 



Figure 2 presents the graph of the Lyapunov exponent logA r as a function of r. 

Lemma 5.2 Let z be a solution of (4-9) such that \e z \ ^ 1 and e z G R. Then dF/dz ^ 
(in fact, dF/dz is a real negative number). 

Proof. A direct calculation shows that 

dF_ _ (l+r + r 2 )(l+e z ) 2 + r 3 (l-e 2 ) 2 
~dz~ ~ 4r 2 (l +r)(l -e z ) 2 

which proves the lemma. □ 

For any r = m/n with an even n and odd m, the corresponding solution e z = — A r 
satisfies the conditions of the above lemma. Hence, this solution changes continuously 
with r, and so we get the following 

Corollary 5.3 For every r = m/n with an even n and an odd m there is an interval 
(r — e, r + e) in which all parameter values have unstable perturbations. 

Therefore, unstable perturbations exist for an open and dense set of parameter values 
< r < 1. One would naturally wonder if all r's have unstable perturbations. The 
answer is, surprisingly, negative: 
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Fact 5.4 For the density (5.1) with r = 1/3, there are no solutions z of (4-9) with Re(z) > 
0, hence there are not solutions of the linearized equation which grow exponentially with 
time. 

Proof. The characteristic equation 



has two real roots (A = ±1) and two complex roots. The complex roots are, on the one 
hand, conjugate and, on the other, satisfy the reciprocability rule: -P(A) = if and only 
if P(l/A) = 0. Hence, they must belong to the unit circle |A| = 1. □ 

It is interesting to know if other rational parameter values r = m/n with odd n are 
also stable. We have examined the values r = 1/n for small odd values of n — 5, 7, . . . , 31 
numerically (by using MATLAB) and always found that all the roots of -P(A) belonged to 
the unit circle. Therefore, we conjecture that the values r = 1/n with odd n are stable. 

On the other hand, the values r = m/n with an odd n but m > 1 appear to be 
unstable. For r = 2/3, 2/5, 3/5, 3/7 we found, again numerically (by using MATLAB), 
roots A such that |A| > 1. All those roots are complex, for example for r = 2/3 they are 
A = —0.3778 ±1.7173 i. It remains to determine theoretically whether all rational values 
r = m/n with m > 1 are unstable, we leave this question open. 

Fact 5.4 seems to disagree with Theorem 5.1. Indeed, let po(v) be the rectangular 
density (5.1) corresponding to r = 1/3 and p(x, v,0) = Po{v) + epi(x,v,0) an arbitrary 
perturbation with an infinitesimally small e. According to Fact 5.4, this perturbation 
cannot grow exponentially in time. On the other hand, let us approximate 1/3 by a 
rational number r = m/n with an even n. Denote by Pq(v) the corresponding rectangular 
density (5.1) for the chosen r = m/n. Then we have 



Hence, if |r — 1/3| < e, then (p —p^/e is of order one (in the L l metric), and p(a;, t>, 0) 
becomes an e-perturbation of the density Pq(v). As such, it "must" grow exponentially 
in time according to Theorem 5.1. This apparent disagreement requires an explanation, 
which we provide next. 

We recall that smooth unstable perturbations are given by the general formula (4.8). 
For the rectangular density (5.1), the velocity v in this formula only takes two values, 
\v\ — r and |i>| = 1, hence the factor C(\v |) takes two values, as well, and so plays little 
role. For simplicity, we set \v | = 1 and ignore the constant factor C(|t>|) = C(l). Now 
the (real part of) unstable perturbations is described by 



A similar formula holds for \v\ = r, and we omit it. Now recall that for any rational r = 
m/n we have e z l m = A, where A < —1 is the eigenvalue of B described by Theorem 5.1. 
Therefore, Rez = mlog|A| = log A and Imz = ±mn. 



p(x,v,0) =p* (v) +ep 2 (x,v,0) with p 2 = pi + (p - p* )/s 



(5.8) 



q(y, 1, t) = Re e z{t - y) = e iRcz)it - y) cos[(Im z) (t - y)\ 



(5.9) 



20 



We see that the real part of z changes continuously with r = m/n but the imaginary 
part does not. In particular, when r = m/n is close to 1/3 and n is even, both m and n 
have to be large, so that m — > oo and | Im z\ — ► oo as r — > 1/3. In terms of the perturba- 
tion (5.9), the growth of |Imz|, as r approaches 1/3, implies that the function q(y, l,t) 
becomes highly oscillatory, and so does the corresponding initial unstable perturbation 
p(x,v,0) = h(v)q(x,v,0). Thus, the linear subspace of unstable perturbations (along 
which exponential growth takes place) becomes nearly orthogonal to any given function, 
in particular to p%{ x i v, 0) defined in (5.8). 

This explains the above "disagreement" . The density P2 does grow exponentially in 
time for any r = m/n with an even n, but, as r — > 1/3, the projection of P2 onto the 
unstable subspace corresponding to the positive Lyapunov exponent logA r > becomes 
small and vanishes in the limit, hence the exponential growth is not visible during a long 
initial interval of time. In the limit r — > 1/3 that "initial interval" becomes infinite and 
the instability evaporates. 

One can also reverse this line of argument. Indeed, when e is not infinitesimally 
small but finite, the representation (5.8) implies that any perturbation p(x, v, 0) of the 
rectangular density (5.1) for any < r < 1 will eventually grow exponentially fast in 
time (because any r e (0, 1) can be approximated by rational numbers m/n with even 
n). We checked this conclusion experimentally and found that it was indeed correct. 



k ^ "bump" 




1 






X 




oTs 1 











Figure 3: Initial rectangular density (5.1) perturbed by a "bump". 

To investigate the instability experimentally, we solved the hydrodynamical equations 
(H1)-(H4) numerically starting with a perturbed rectangular density (5.1) shown on 
Fig. 3. The initial density p(x, v, 0) takes the value one on the black region and zero 
elsewhere. The small 'bump" on the top left edge of the upper rectangle represents the 
perturbation. The area of the bump in our experiments was less than 1% relative to the 
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total area of each black rectangle. 
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Figure 4: Piston's trajectory for a perturbed rectangular density with r = 1/3. 
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Figure 5: Piston's trajectory for a perturbed rectangular density with r = 1/4. 



Figures 4 and 5 show typical trajectories of the piston X(t) for r = 1/3 and r = 1/4, 
respectively. One can see that by the time t = 15 the piston's motion shows signs of 
exponential instability in both cases, but there is a notable difference. For r = 1/4 the 
piston just swings back and forth with a monotonically increasing amplitude, as time goes 
on. For r = 1/3, the amplitude of oscillations grows slower but the frequency increases 
quickly. The higher frequency of oscillations of X(t) for r = 1/3 probably reflects the 
oscillatory structure of unstable perturbations for rational r = m/n approximating 1/3. 



6 Periodic solutions of the hydro dynamical equations 

Here we discuss the long-term behavior of our system in the unstable regime. 
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In our previous work [CL] we reported the results of computer simulations of the 
piston dynamics in an ideal gas with many (up to 27 million) particles. The initial 
configuration of particles was selected randomly with the average density (4.7), see [CL] 
for details. A typical trajectory of the piston X(t) found in our experiments is shown 
here on Fig. 6. One can see that during the initial interval of time < t < 8 the piston 
moves back and forth with an exponentially increasing amplitude, which is consistent 
with our analysis in Section 4, where the density (4.7) was proven to be unstable. 




5 10 15 20 25 30 



Figure 6: Piston's trajectory in the mechanical model with 10 6 particles. 

Later on, however, at times 8 < t < 15, the amplitude of the piston's oscillations 
decreases to a certain constant value (nearly a half of its maximum attained at t — 8). 
Then the piston's oscillations become very stable and continue almost unchanged for a 
very long time, up to t = 50 or 100, with a very slowly decreasing amplitude. 

On the other hand, we have solved the hydrodynamical equations (H1)-(H4) numer- 
ically, starting with the same initial density (4.7) perturbed by a bump shown on Fig. 3. 
Figure 7 presents the resulting trajectory of the piston. One can see that it behaves al- 
most identically to the simulated trajectory of the piston shown on Fig. 6. Thus, not only 
the initial instability, but also the long term behavior of the simulated piston trajectory 
match those of perturbed solutions of the hydrodynamical equations. 




Figure 7: Piston's trajectory from the HE for a perturbed rectangular density with 

r = 1/2. 
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The behavior shown on Fig. 7 persists when various perturbations of the initial density 
(4.7) are applied. It seems like there is a periodic cycle or an invariant manifold of 
quasi-periodic solutions of (H1)-(H4) that acts as an attractor. Of course, due to the 
time-reversibility of the hydrodynamical equations there can be no attractors in the strict 
sense. It is more likely that there is an invariant manifold of periodic or quasi-periodic 
solutions that acts as a saddle point in the phase space: typical trajectories approach 
that manifold temporarily and then slowly move away. We cannot rigorously prove the 
existence of periodic or quasi-periodic solutions, but we construct such solutions by using 
perturbative analysis. 

We will be looking for solution of the hydrodynamical equations (H1)-(H4) such that 
the piston makes harmonic oscillations 

X(t) = - + ecosut, X(t) = — euo sincut, (6.10) 

with some fixed u> > and small e > 0. We will approximate such solutions up to the 
first order in e, i.e. ignoring terms of higher order. 

The construction is done in two steps. First, we assume that the piston moves as 
prescribed by (6.10) and consider the motion of a fluid point bouncing against the moving 
piston X(t) and the fixed wall x — 0. Second, we define the density p(x,v,t) which, 
coupled with the piston's oscillations (6.10), satisfies equations (H1)-(H4). 

Let the piston move according to Eqs. (6.10). Then fluid points in the left compart- 
ment < x < X(t) bounce against the wall x = and the piston, the latter simply acts 
on them as a moving wall. It is known that the phase space of gas particles bouncing 
against a periodically moving wall necessarily contains many invariant curves. Moreover, 
the region corresponding to high velocities \v\ > Vq is densely filled by such invariant 
curves, the larger vq the higher the density of invariant curves. This fact is a conse- 
quence of KAM theory, it was first proved by R. Douady in his thesis [Do] and later 
independently by S. Laederich and M. Levi [LL]. We describe these invariant curves 
approximately, up to the first order in e, by equation 

v + eF(t, V) = V + 0(e 2 ) (6.11) 

where v denotes the velocity of the particle when it kicks the piston, t the the collision 
time, and V is the parameter of the curve. In fact, we will construct invariant curves for 
all V > Vq with some Vq > 0. Here we only consider particles to left to the piston, the 
particles to the right of the piston are completely symmetric. 

Let us consider successive collisions of a gas particle with the piston. Denote by 
v n > the velocity of the particle before its n-th collision and by t n the time of that 
collision. Then the law of elastic impact reads 

v n+ i = v n - 2X(t n ) =v n + 2eusmut n (6.12) 
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Let t n+ i/2 denote the time at which the particle bounces off the wall x = between its 
n-th and (n + l)-st collisions with the piston. Obviously, t n +i/2 = t n + X(t n )/v n+ i and 
v n+ i(t n+ i — t n+ i/2) = X(t n+ i). Since we are interested in knowing v n up to terms 0(e), 
it is sufficient to find t n up to terms C(l). This is easy: 

t n +i = t n + — + 0(e) = t n + i + 0(e) (6.13) 

V n +1 V 

where we used (6.11). 

Now let us look for an invariant curve of the form 

v + eF(t,V) = V + 0(e 2 ) 

We have to impose the constraint 

v n+l + eF(t n+1 , V)=v n + eF{t n , V) + 0{e 2 ) (6.14) 

By equations (6.12), (6.13) and (6.14) we get 

v n + 2etusmujt n + eF{t n + 1/V + 0(e), V) = v n + eF(t n , V) 

Cancelling v n and e and removing the index n gives a general equation for an invariant 
curve: 

2ujsmut + F(t + 1/V,V) - F(t,V) = (6.15) 
We construct solutions of this equation in the form 

F(t, V) = a cos ut + bsinut (6.16) 

where a and b depend on V. By substituting this expression into (6.15) we find that 
(6.15) can only hold if 

a(cos(u/V) - 1) + bsm(u;/V) = 
asm(uj/V) — b(cos(uj /V) — 1) = 2u 

The solution of the above system is 

oj sin(cj/V) 



1 — cos(c<j/V) 
b = oj 

Remark. Notice that a (and hence the invariant curve) is not defined for V = 

k = ±1, ±2, ... To avoid these singularities, we will not use invariant curves corresponding 

to V < j-. In particular, the density p(x, v,t) that we define below will be constant for 
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Thus, for any V > Vo > ^ we can define an invariant curve u(x, t] V) in the phase 
space of gas particles, where V is the parameter of the curve and u(x, t; V) is the velocity 
of the particle on the curve at point x at time t. The curve is made by two branches: 
the upper branch u + and the lower branch u~ . Obviously, we have 

u + (X(t),t;V) = V-eF(t,V) + 0(e 2 ) 
u-(X(t),t;V) = -[V -eF(t,V)-2X(t)] + 0(e 2 ) 
u + (t,0,V) = -u-(t,0,V) 

Note that the last equation here is equivalent to (6.14). 

Now we define a density p(x, v, t) so that its value on each invariant curve u(x, t; V), 
\V\ > V , is a constant denoted by p(V). Between the curves u + (x,t; Vq) and u~(x,t; V ) 
we set the density to a constant equal to one. Therefore 

p(x,u + (t,x,V),t) =p(x,u~(t,x, V),t) = p(V) if V >V 

and 

p(x,v,t) = l if u~(x, t; Vq) < v < u + (x, t; V ) 
The function p(V) and the "cutoff" value V > ^ will be specified below. 
Example. Let us set p(V) = for V > V , i.e. 



p(x,v,t) 



1 for u (x, t] Vo) < u < u + (x, t] Vq) 
elsewhere 



In order to compute the pressure on the piston we only need to know the density p(x, v, t) 
at the point x — X(t), i.e. we need to know the function 

v(t, V) := u + (X(t),t, V) = V + eF(t, V) + 0{e 2 ) 

In our example the density on the piston (on the left hand side) is 1 up to v + = Vq — 
eF(Vo, t). The density on the piston on the right hand side is 1 up to a similar invariant 
curve, which is phase shifted by At = n/uo. Therefore the density on the right hand side 
is 1 down to v~ = -V - eF(t, V ). Since F(t + tt/uj) = -F(t) by (6.16), the velocity of 
the piston is exactly the average of v + and v~ and therefore is —eF(t, Vo). 

Thus our density and the piston satisfy the hydrodynamical equations (H1)-(H4) if 
X = —eF(t, Vo), which gives 

eoj sin(a;/Vo) cosut 

—euj sin out = —eu sin cut — : — - — 

1 — cos(cj/Vo) 

In our example, the only possible choice is Vo = ^. 

Now let us consider the case of a generic function p(V). The pressure on the piston 
on the left hand side is equal to 

Pl= / Pl(v) (v - Xf dv = / p L {v){v 2 -2vX)dv + 0{e 2 ) 
Jx Jo 
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where pl(v) = p(X(t) — 0,v,t) is the density on the piston (we have used the fact that 
X = 0(e)). Recall that the density is p L (v) = p(V) = p(v + eF(t,v)) + 0{e 2 ). /From 
now on we neglect terms of order 0(s 2 ). Then we get 

roo 

P L = (v 2 -2vX) p(v + sF(t,v))dv 
Jo 

The pressure on the right hand side is given, by analogy, 

/0 roo 
p R (v) {v 2 - 2vX) dv = / p R (-v) {v 2 + 2vX) dv 
-oo JO 

Note that for v > we have pr{v) = p(V) = p{v + eF(v,t + 7r)) = p{y — eF(t,v)). 

Therefore _ 

P R = (v 2 + 2vX)p(v -eF(t,v))dv 
Jo 

We now conclude that Pl = Pr iff 

■ = fS°pf(v)F(t,vWdv _ fS°(/(v)F(t,v)v*dv 
£ J °°p(v)2vdv ~ 6 f °°p>(v)v 2 dv 

which is analogous to our early formula (4.1). 

Using (6.16) and the subsequent equations we find 

roo // x sjMw/v) 2 J 

. , Jo Pv v ) i-cos^ao v av , 

X = —euj sin out — e — uj cos uot 

J p'(v)v 2 dv 

Our density, coupled with the piston oscillations (6.10), satisfies the hydrodynamical 
equations (H1)-(H4) if and only if X = —ecu smut. This implies 



J<j, 



dvp\v)v 2 S[niu ! V ] , =0 (6.17) 

where we have imposed p' = for v < u/2n. 

Interestingly, (6.17) seems to be related to our early equation (4.9). Precisely, let 
z in (4.9) be a purely imaginary number, z = ivi. Also note that p(v) in (4.9) is just 
proportional to v 2 p'(v) here. Then (6.17) becomes equivalent to ImF(z) = 0, with F(z) 
defined by (4.10). In other words, (6.17) expresses the "imaginary part" of the equation 

(4.9) . We already observed in the previous section that Imz characterized the frequency 
of oscillations of unstable perturbations, and here cu = Im z is the frequency of oscillations 
of the piston. We note that for z = ui one always has ReF(z) = 0, as it follows from 

(4.10) , hence in our case (6.17) is equivalent to F(z) = 0. 

Next, we note that p'(y) must be supported on both the intervals (u/2n, uj/tt] and 
[cu/ir, +oo). In fact the fraction in Eq. (6.17) is negative for v G (u/2tt, uj/tt) and positive 
for v > u/n. 
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Obviously Eq. (6.17) can be satisfied in many different ways, all of them leading to 
different solutions of the hydrodynamical equations. 

Lastly, we want to examine a special case when the initial density p(x, v, 0) = p(y) 
is a monotonic function in \v\, i.e. when the hydrodynamic equations are stable. Then 
there are some quantitative restrictions on the period of oscillations of the piston. The 
period T = 2n/u can be bounded from below by a function of the average kinetic energy 
(K) = K/M, where 

poo ijj^ roo 

K— p( v )tt^ v an d M— p(v) dv 
Jo 2 Jo 

Proposition 6.1 Let p' < be supported on the interval [uj/2h, oo) and satisfy Eq. 
(6.17). Then the period of oscillations T is bounded by 



T > 



3(K) 



The equality holds when p' = {n/u))5{v — u/n), i.e. when p{\v\) is constant on (0,u/ir) 
and elsewhere. 

Even though the proposition is stated for monotonic densities only, let us apply it 
to our unstable density (4.7). In this case K = 7/24 and T = ^ff ~ 1.51186. The 
experimentally determined period of oscillations of the piston is T ~ 1.62, see [CL]. We 
also simulated the piston trajectory with other unstable densities (5.1) with r — > and 
observed that the period of oscillations approached its lower bound 2 given by the above 
proposition. 

Proof of Proposition 6.1. Consider a function 

sm(u/v) 

G(v) := v — — 

1 — cos{u/v) 

Then equation (6.17) reads 

C := / -p'vGdv = 

Juj/2-k 

Note that in the interval [u/27r, oo) the function G(v) is strictly increasing and that 
G{u)/2n) = — oo, G(u/n) = 0, and G(oo) = oo. 

Introducing a new function R(v) = —p'(v)v > and integrating by parts yields 

M= / R(v)dv, K= / R(v)—dv, C= / R(v)G(v) dv 

Jw/2-k Jw/2-k 6 Jw/2-k 

It is useful to replace v by a new variable u = G(v), — oo < u < oo. Since G is strictly 
increasing, we can write 

/oo roo roo 

S(u)du, K= / S(u)r)(u)du, C= / S{u)udu 
-oo J —oo J oo 
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where S(u) = RiG' 1 (u)) / G' (G' 1 (u)) and rj(u) = (G-\u)) 2 /6 (here G~ l denotes the 
inverse of the function G). 

We have to solve the following variational problem: minimize K/M under the con- 
straint C = 0. As we shall see in the sequel the function r\ turns out to be convex. This 
easily implies Proposition 6.1. 

Also, the convexity of i] implies that the solution of the variational problem S is a 
delta-function centered at u — 0, i.e. at v = uj/tt (so that C vanished). 

So it only remains to prove that i] is a convex function of u. By direct computation 
we get 

Q V ' = ^(G~\u)f = 2G-\u) dG ^ u) = 2v/G'(v) 
du du 

Hence it is sufficient to prove that the function G'{v)/v is strictly decreasing in the 

interval v > uo/2n. Without loss of generality we set u — 1, then 



1 — cos(l/f ) 
Consider a new function 

G'(v) l + usin(l/u) 



H(v) 



v v 2 (l — cos(l/f )) 
then 

, y — v — (— 1 + v 2 ) sm(l/v ) + v cos(l/w)(l + v sin(l/f )) 
= t; 4 (l - cos(l/t;)) 2 

The denominator of H' being positive, we only need to prove that the numerator of H' 
is negative in the interval v > l/2n. 

If we replace v by 1/x and multiply by the numerator by x 2 we find the expression 

h(x) = — x + (x 2 — 1) sin(x) + cos(a;)(a; + sinx) 
= (cosx — l)(sina; + x) + x 2 sin a; 

We need to show that h(x) < in the interval x E (0, 2ir). First of all, h(0) = h(2n) = 
and for any x G (it, 27r) the expression is clearly negative. 

It only remains to prove that h(x) is negative in (0, 7r]. By computing the Taylor 
expansion of h about h — one finds 

+00 r)2fc _ A h,2 

It is easy to prove that for any x G (0, it] this is an alternating series, the absolute values 
of its terms being strictly decreasing. 

The first few terms of the above expansion are 

h(x) = -— (l--x 2 + —x A - -^—x 6 ) + 0(x 15 ) 
y ' 180 V 21 240 166320 / v ' 
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Therefore 

x 7 . x 7 ( 2 2 \ 

< h(x) < 1 x 2 

180 v ; 180 V 21 J 

which implies that h(x) < for any x G (0, rr]. □ 
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